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ABSTRACT.  We  discuss  an  adaptive  local  refinement  finite  element  . 
method  for  solving  initial-boundary  value  problems  for  vector  systems  of 
partial  differential  equations  in  one  space  dimension  and  time.  The  method 
uses  piecewise  bilinear  rectangular  space-time  finite  elements.  For  each  time 
step,  grids  are  automatically  added  to  regions  where  the  local  discretization 
error  is  estimated  as  being  larger  than  a  prescribed  tolerance.  We  discuss 
several  aspects  of  our  algorithm,  including  the  tree  structure  that  is  used  to 
represent  the  finite  element  solution  and  grids,  an  error  estimation  technique, 
and  initial  and  boundary  conditions  at  coarse-fine  mesh  interfaces.  We  also 
present  computational  results  for  a  simple  linear  hyperbolic  problem,  a 
problem  involving  Burgers'  equation,  and  a  model  combustion  problem. 


1.  INTRODUCTION.  There  is  an  ever  increasing  need  to  solve 
problems  of  greater  complexity  and  a  corresponding  need  for  reliable  and 
robust  software  tools  to  accurately  and  efficiently  describe  the  phenomena. 
Adaptive  techniques  are  good  candidates  for  providing  the  computational 
methods  and  codes  necessary  to  solve  some  of  these  difficult  problems.  Two 
popular  adaptive  techniques  are:  (i)  moving  mesh  methods,  where  a  grid  of  a 
fixed  number  of  finite  difference  cells  or  finite  elements  is  moved  so  as  to 
follow  and  resolve  local  nonuniformities  in  the  solution,  and  (ii)  local 
refinement  methods,  where  uniform  fine  grids  are  added  to  coarser  grids  in 
regions  where  the  solution  is  not  adequately  resolved.  A  representative 
sample  of  both  types  of  methods  is  contained  in  Babuska,  Chandra,  and 
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Flaherty  [2] .  Recently,  Adjerid  and  Flaherty  [1]  deveiopecf* element 
method  that  combines  mesh  moving  and  refinement. 

Herein,  we  discuss  a  local  refinement  finite  element  procedure  for 
finding  numerical  solutions  of  M-dimensional  vector  systems  of  partial 
differential  equations  having  the  form 


Lu  :=  ut  ♦  f(x,t,u,ux)  -  lD(x,t,u)ux]x  *  0, 


a  <  x  <  b,  t  >  0, 


subject  to  the  initial  conditions 
u(x,0)  =  u*Cx)  ,  a  £  x  £  b  , 


(1.1) 


(1.2) 


and  appropriate  boundary  conditions  so  that  the  problem  has  a  well  posed 
solution. 

We  discretize  (1.1,2)  for  a  time  step  using  a  finite  element-Galerkin 
procedure  with  piecewise  bilinear  approximations  on  a  rectangular  space-time 
net.  At  the  end  of  each  time  step  we  estimate  the  local  discretization  error, 
add  finer  subgrids  of  space-time  elements  in  regions  of  high  error,  and 
recursively  solve  the  problem  again  in  these  regions.  The  process  terminates 
when  the  error  estimate  on  each  grid  is  less  than  a  prescribed  tolerance.  The 
original  coarse  space-time  grid  is  then  carried  forward  for  the  next  time  step 
and  the  strategy  is  repeated.  Our  algorithm  is  discussed  further  in  Flaherty 
and  Moore  [9]  and  some  of  this  discussion  is  repeated  in  Section  2. 

Berger  [3]  used  a  similar  local  refinement  procedure  to  solve  one-and 
two-dimensional  hyperbolic  systems.  She  used  explicit  finite  difference 
schemes  to  discretize  the  partial  differential  equations,  while  we  use  implicit 
finite  element  techniques  since  we  are  primarily  interested  in  parabolic 
problems. 

In  addition  to  the  discretization  technique,  the  major  numerical 
questions  that  must  be  answered  as  part  of  the  development  of  a  local 
refinement  code  are  (i)  the  estimation  of  the  discretization  error  and  (ii)  the 
appropriate  initial  and  boundary  conditions  to  apply  at  coarse-fine  mesh 
interfaces.  Of  course.  Computer  Science  questions,  such  as  which  language 
to  use  to  describe  and  implement  the  various  algorithms  and  what  data 
structures  to  use  to  represent  and  store  the  grids  and  solutions  must  also  be 
answered.  Our  work  in  all  of  these  areas  is  still  far  from  complete  and  herein 
we  only  discuss  our  progress  and  thoughts  on  error  estimation  techniques, 
data  structures,  and  interface  conditions  (cf.  Section  2).  In  Section  3,  we 
present  the  results  of  three  examples  that  illustrate  our  method  and  the 
discussion  of  Section  2,  and  in  Section  4,  we  present  some  preliminary 
conclusions  and  future  plans. 


2.  FINITE  ELEMENT  ALGORITHM.  We  discretize  equation  (1.1)  on  a 
strip  a  <  x  <  fl,  p  <  t  <  q  using  a  finite  element-Galerkin  method  with  a 
uniform  grid  of  N  rectangular  elements  of  size  (6  -  o)/N  by  (q  -  p).  We 
refer  to  this  grid  as  R(a.3,p.q,N,f,s),  where  f  and  s  are  pointers  to  the 
father  and  son  grids  discussed  later.  Each  grid  uses  records  to  store  the 
appropriate  information. 


We  generate  the  discrete  system  on  R(a,0,p,q,N,f,s)  in  the  usual 
manner;  thus,  we  approximate  u  by  U(x,t)  and  select  test  functions  V(x,t), 
where  U  and  V  are  elements  of  a  space  of  C*  bilinear  polynomials  with  respect 
to  the  grid  R.  We  then  take  the  inner  product  of  equation  (1.1)  and  V, 
replace  u  by  U,  and  integrate  any  diffusive  terms  by  parts  to  obtain 

/  (VTUt  ♦  VTf(x,t,U,Ux)  ♦  V^D(x,t,U)Ux]dxdt 


q  t 

-/  V  D(x,t,U)U  |pdt  =  0  .  (2.2) 

x  s 

P 

Equation  (2.2)  must  vanish  for  all  bilinear  functions  V  on  the  grid  R.  The 

integrals  are  approximated  using  a  four  point  Gauss  quadrature  rule  and  the 

resulting  nonlinear  system  is  solved  by  Newton  iteration  (cf.,e.g.,  [7]  for 

additional  details).  Appropriate  initial  and  boundary  conditions  for  (2.2)  are 
discussed  later  in  this  section. 

We  describe  our  local  refinement  procedure  for  solving  problem  (1.1,2) 
for  one  time  step  (t*,tl)  on  a  coarse  grid  with  N*  elements,  i.e,  on 

R(a,b,t',tl,N*,0,s)  (where  the  pointer  f  =  0  signifies  that  this  grid  has  no 
father).  To  solve  this  problem  we  simply  call  the  procedure  "locref"  with  the 
arguments  Rta.b.tVtVN'.O.s),  tol,  tsub  for  each  coarse  grid  time  interval. 
A  pseudo-PASCAL  description  of  the  procedure  "locref"  is  shown  in  Figure  1. 


procedure  locref  (R(a,8.p,q.N,f,s),  tol,  tsub) 
begin 

Solve  the  finite  element  equations  (2.2)  on  R(a,B,p,q,N,f,s); 

Estimate  the  error  on  R(tx,&,p,q,N,f,s); 
if  error  >  tol  then 
begin 

calculate  where  error  >  tol  and  return  the  son  grids; 
for  j  :*  1  to  tsub  do 

for  i  :«  1  to  number  of  sons  do 
begin 

P(j]  :=  P  *  (j-1)*(q-p)/tsub; 
q[jj  :*  Ptj]  *  (q-p)/tsub; 
locref  (R(a[i],6[i],p[j],q[j],N[i], 
R(a,S,p,q,N,f,s),s[iJ,  tol,  tsub) 

end 

end 

end; 

Figure  1.  Algorithm  for  local  refinement  solution  of  (1.1,2)  on 
R(a,8,p,q,N,f,s)  with  an  error  tolerance  of  tol  and  dividing  the  local  time 
step  by  tsub  each  time  the  error  test  is  not  satisfied. 


The  recursive  algorithm  locref  sets  up  a  tree  structure  of  grids  with 
R  •:  a  .  o .  t5 ,  t ;  N  * ,  0.  s)  being  the  root  node  and  with  the  solution  being 


generated  by  a  preorder  traversal  of  the  tree  at  each  local  time  step.  For 
example,  if  the  root  grid  is  refined  to  give  two  subgrids  and  the  time  step  is 
halved,  then  the  problem  is  solved  on  the  first  subgrid  on  its  first  time  step, 
then  on  the  second  subgrid  on  the  same  time  step,  then  this  procedure  is 
repeated  for  the  second  time  step.  The  error  is  estimated  by  Richardson 
extrapolation,  i.e.,  the  space  and  time  steps  are  halved  and  the  problem  is 
solved  again  on  this  new  grid.  The  two  solutions  that  are  obtained  at  each 
original  grid  point  are  used  to  generate  an  error  estimate.  If  this  pointwise 
estimate  exceeds  the  tolerance  "tol”,  finer  grids  are  added  as  leaf  nodes  to 
the  tree.  This  procedure  is  similar  to  one  used  by  Berger  [3);  however, 
there  are  more  economical  error  estimation  strategies  (cf.,  e.g.,  Bieterman 
and  Babuska  [5,  6])  which  we  are  currently  investigating. 

In  order  to  solve'  the  finite  element  system  (2.2)  we  need  to  supply 
initial  and  boundary  conditions.  On  any  grid  with  p  -  0,  e  =  a,  or  ft  =  b 
these  can  be  obtained  from  the  initial  condition  (1.2)  or  prescribed  boundary 
conditions.  However,  artificial  initial  and  boundary  conditions  must  be  created 
at  all  other  coarse-fine  mesh  interfaces.  This  is  a  difficult  and  crucial 
problem  that  is  discussed  for  explicit  finite  difference  methods  by  Berger 
[3,  4];  however,  it  is  largely  unanswered  for  finite  element  applications. 
Instabilities  or  incorrect  solutions  (cf.  Example  1  of  Section  3)  can  result  if 
inappropriate  conditions  are  specified. 

For  initial  conditions,  two  strategies  immediately  come  to  mind:  (i) 
saving  all  fine  grid  data  for  propagation  in  time  or  (ii)  interpolating  the  best 
coarse  grid  data  to  finer  grids.  We  consider  a  blend  of  the  two  strategies 
which  consists  of  saving  the  fine  grid  data  down  to  a  given  level  X  in  the 
tree  and  subsequently  interpolating  for  finer  grids.  Each  grid  in  the  first  X 
levels  either  has  a  linked  list  of  the  initial  data  directly  associated  with  it  or 
uses  an  initial  data  list  of  an  ancestor  grid.  To  find  the  value  of  the  solution 
at  some  new  initial  point,  the  coordinate  of  that  point  is  sequentially  compared 
to  values  in  the  linked  list  until  an  interval  containing  the  point  is  found  so 
that  interpolation  can  be  used.  This  is  costly  and  we  are  investigating  more 
efficient  procedures  that  use  the  natural  ordering  that  already  exists.  We 
used  either  piecewise  linear  interpolation  or  piecewise  parabolic  interpolation 

with  shape  preserving  splines  developed  by  McLaughlin  [10].  For  each  grid 

in  the  first  X  levels  of  the  tree,  a  linked  list  is  created  to  store  the  initial 
data.  We  are  studying  several  alternative  ways  of  determining  a  proper  value 
for  X. 

At  the  present  time,  we  prescribe  internal  Dirichiet  boundary  conditions 
by  linearly  interpolating  from  coarse  to  finer  grids.  A  buffer  zone  of  two 

elements  is  added  to  each  end  of  regions  of  high  error  that  do  not  intersect 

the  boundaries  x  *  a  and  b.  If  two  buffer  zones  overlap  or  are  separated 
from  one  another  by  one  element,  the  two  grids  are  joined.  Similarly,  if  the 
buffer  is  only  one  element  away  from  either  a  or  b,  that  element  is  added  to 
the  grid. 

3.  NUMERICAL  EXAMPLES.  An  experimental  code  based  on  the 
algorithms  in  Section  2  has  been  written  in  FORTRAN-77.  We  are  testing  it 
on  several  examples,  some  of  these  follow  and  others  are  presented  in  [9]. 
All  results  were  computed  in  double  precision  on  an  IBM  30810  computer. 


Example  1.  In  order  to  illustrate  the  importance  of  adequately  resolving 
initial  conditions  at  each  time  step  we  solve  the  linear  hyperbolic  initial  value 
problem 


(( 1/2)  ( cos (20ir(x -0.45))  -  1)  , 

u(x,0)  =  u* (x)  =  j  0.35  <  x  <  0.75 

10  ,  otherwise 

We  solve  this  problem  for  one  coarse  time  step  of  At  -  0.05,  10  elements  on  0 
<  x  <  1,  tol  =  0.01.  For  small  enough  times  the  exact  solution  is  u*(x-t).  If 
initial  conditions  are  interpolated  from  the  coarse  to  the  fine  grid,  the 
oscillations  are  missed  and  ao  incorrect  solution  is  computed,  possibly  without 
a  user  realizing  that  there  is  anything  wrong.  However,  saving  initial  values 
for  the  first  8  levels  of  the  tree  of  grids  calculates  the  correct  solution  to  the 
prescribed  accuracy.  The  incorrect  -and  correct  solutions  are  shown  at 
t  2  0.05  in  Figure  2. 

Example  2.  We  solve  the  following  problem  for  Burgers'  equation: 
ut  •  uux  =  du>x  ,  0  <  x  <  1  ,  0  <  t  <  t  , 

u(x,0)  =  sinirx  ,  0  <  x  <  1  , 

u(0,t)  =  u(1,t)  2  0  ,  t  >  0  . 

We  choose  d  =  0.00003,  a  coarse  grid  of  10  elements  and  At  s  0.1,  and 
piecewise  parabolic  approximations  for  the  initial  conditions  with  X  *  6.  It  is 
well  known,  that  the,  solution  of  this  problem  is  a  "pulse”  that  steepens  as  it 
travels  to  the  right  until  it  forms  a  shock  layer  at  x  2  1.  After  a  time  of 
0(l/d)  the  pulse  dissipates  and  the  solution-  decays  to  zero.  We  solve  this 
problem  for  tol  *  0.01  and  0.001  and  show  the  solutions  at  t  =  0.4  in  Figure 

3.  The  solution  with  the  cruder  tolerance  is  exhibiting  some  oscillations  that 

are  within  our  bounds.  These*„.howev*r,  ~ere~~not  visible  when  the  finer 
tolerance  t»-  used  to  solve  the  problem. 

Example  3.  We  solve  the  model  combustion  problem 
“t  *  u„  •  2.“  *  uxx  .  0  <  x  *  I  .  0  ‘  t  <  1  , 
u(x,0)  =  0  ,  u(0,t)=0  ,  ux(1,t)  *  0  . 

The  exponential  nonlinearity  is  typical  io  combustion  problems  having 
Arrhenius  chemical  kinetics.  However,  in  this  case  the  solution  develops  a 
"hot  spot"  at  x  s  1  and  becomes  infinite  when  t  is  approximately  0.85.  We 
choose  a  coarse  grid  of  20  elements  and  At  s  0.05,  tol  2  0.001,  and  X  2  6.  In 
Figure  4  we  show  the  computed  solution  U(x,t)  as  a  function  of  x  for  t  2 
0.05,  0.6,  and  0.8  and  in  Figure  5  we  show  the  mesh  that  was  used  to  solve 
the  problem.  We  see  that  the  mesh  is  initially  concentrated  in  the  region  near 
x  2  0  where  the  curvature  of  the  solution  is  largest.  As  time  progresses  and 
the  curvature  diminishes,  excessive  refinement  is  not  necessary.  Finally,  as 
the  solution  begins  to  "blow-up"  our  algorithm  generates  a  fine  mesh  only  in 
the  region  near  x  =  1 . 


4.  DISCUSSION  AND  CONCLUSIONS.  We  have  briefly  described  an 
adaptive  local  refinement  algorithm  for  solving  time  dependent  partial 
differential  equations.  Even  though  this  is  very  much  a  working  algorithm, 
and  not  a  production  code,  we  are  very  encouraged  by  the  preliminary 
results.  We  are  investigating  several  possible  ways  of  improving  the 
efficiency  and  robustness  of  our  algorithm.  These  include  adding  higher 
order  '  polynomial  finite  element  approximations,  adaptively  changing  the 
number  of  elements  that  are  carried  forward  in  the  coarse  grid  at  each  coarse 
time  step,  how  to  select  the  appropriate  buffer  length,  adaptively  determining 
the  optimal  number  of  levels  of  initital  conditions  to  keep  at  coarse-fine 
interfaces,  and  the  best  boundary  conditions  to  apply  at  internal  boundaries. 
We  are  encouraged  by  the  performance  of  McLaughlin's  [10]  shape  preserving 
parabolic  splines;  however,  the  entire  area  of  interpolating  from  coarse  to 
fine  grids  needs  further  study.  We  are  also  developing  non-Dirichlet 
"natural"  boundary  conditions  to  use  at  coarse-fine  mesh  interfaces. 

Finally,  we  are  very  interested  in  combining  the  moving  mesh  strategy 
of,  e.g.,  [7,  8]  with  the  present  local  refinement  strategy  and  extending  our 
methods  to  two  and  three  dimensions. 
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Figure  2.  Solution  of  Example  1  at  time  t  =  0.05  using  interpolation  from 
the  coarse  grid  to  the  fine  grid  (top)  and  saving  the  initial  values  for  the 
first  8  levels  of  the  tree  (bottom).  The  upper  solution  overlooks  the 
oscillations  and  is  incorrect. 
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Figure  5.  The  grids  generated  in  solving  Example  3  for  0  <  t  <  0.8.  The 
initial  coarse  mesh  has  20  elements  with  At  =  0.05. 


